%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
params = {"ytick.color" : "c",
"xtick.color" : "c",
"axes.labelcolor" : "c",
"axes.edgecolor" : "c",
"text.color" : "c"}
plt.rcParams.update(params)
import random
import math
from copy import deepcopy
import cv2
#Load the image
plt.figure(figsize = (10,10))
img = cv2.imread('lena_gray.jpg',cv2.IMREAD_GRAYSCALE)
plt.imshow(img, cmap='gray')
def snr(oryg, mod):
oryg = oryg.astype('int64')
mod = mod.astype('int64')
a = np.sum(np.square(mod))
b = np.sum(np.square(np.subtract(mod, oryg)))
return a/b
def psnr(oryg, mod):
oryg = oryg.astype('int64')
mod = mod.astype('int64')
a = np.max(mod)
b = np.sum(np.square(np.subtract(mod, oryg)))/(oryg.shape[0]*oryg.shape[1])
return (a*a)/b
def rmse(oryg, mod):
oryg = oryg.astype('int64')
mod = mod.astype('int64')
return np.sqrt(np.sum(np.square(np.subtract(mod, oryg)))/(oryg.shape[0]*oryg.shape[1]))
def add_gaussian(img, sigma):
gaussian_noise = np.random.normal(0, sigma, img.shape)
gaussian_noise = gaussian_noise.reshape(img.shape)
noisy_image = img + gaussian_noise
noisy_image = np.clip(noisy_image, 0, 255)
noisy_image = noisy_image.astype(np.uint8)
return noisy_image
def add_salt_pepper(img, perc):
x, y = img.shape
img2 = deepcopy(img)
for i in range(0,x):
for j in range(0,y):
val = random.random()
if val < (perc/100):
img2[i,j] = random.choice((np.uint8(0), np.uint8(255)))
img2 = np.clip(img2, 0, 255)
img2 = img2.astype(np.uint8)
return img2
def decybel(val):
return 10*math.log(val, 10)
def add_noise(fun, val):
img2 = fun(img, val)
fix, ax = plt.subplots(1,2)
fix.set_size_inches(20,10)
ax[0].imshow(img, cmap='gray')
ax[1].imshow(img2, cmap='gray')
print("RMSE: ", rmse(img, img2))
print("SNR: ",decybel(snr(img, img2)))
print("PSNR: ",decybel(psnr(img, img2)))
return img2
img_salt_pepper_1 = add_noise(add_salt_pepper, 1)
img_salt_pepper_5 = add_noise(add_salt_pepper, 5)
img_salt_pepper_20 = add_noise(add_salt_pepper, 20)
img_salt_pepper_40 = add_noise(add_salt_pepper, 40)
img_gauss_5 = add_noise(add_gaussian,5)
img_gauss_10 = add_noise(add_gaussian,10)
img_gauss_20 = add_noise(add_gaussian, 20)
img_gauss_50 = add_noise(add_gaussian, 50)
def convolution(img, kernel, k, x, y):
value = 0
for i in range(-k, k+1):
for j in range(-k, k+1):
value += kernel[i+k, j+k] * img[x+i, y+j]
return value
def get_area(img, size, x, y):
offset =(size- 1) // 2
start_row = x - offset
start_column = y - offset
area = deepcopy(img[start_row:start_row+size, start_column:start_column+size])
return area
def filter_avg(img, area, center):
x,y = img.shape
img2 = deepcopy(img)
kernel = np.ones((area, area), np.float32)
kernel[area//2, area//2] = center
kernel_sum = np.sum(kernel)
offset =(area- 1) // 2
for i in range(offset, x - offset):
for j in range(offset, y - offset):
img2[i,j] = convolution(img, kernel, offset, i,j) // kernel_sum
return img2
def filter_gauss(img, area, sigma):
kernel = np.zeros((area, area), np.float32)
k = area //2
for i in range(-k, k+1):
for j in range(-k, k+1):
kernel[i+k, j+k] = (1/2*np.pi*sigma**2)*np.exp(-(i**2+j**2)/(2*sigma**2))
kernel_sum = np.sum(kernel)
x,y = img.shape
img2 = deepcopy(img)
offset =(area- 1) // 2
for i in range(offset, x - offset):
for j in range(offset, y - offset):
img2[i,j] = convolution(img, kernel, k, i,j)//kernel_sum
return img2
def filter_median(img, area,val):
x,y = img.shape
img2 = deepcopy(img)
offset =(area- 1) // 2
for i in range(offset, x - offset):
for j in range(offset, y - offset):
img2[i,j] = int(np.median(get_area(img, area, i, j)))
return img2
def filter_conservative_smoothing(img, area,val):
x,y = img.shape
img2 = deepcopy(img)
offset =(area- 1) // 2
for i in range(offset, x - offset):
for j in range(offset, y - offset):
neigh = get_area(img, area, i, j)
neigh[area//2, area//2] = neigh[0,0]
img_max = int(np.max(neigh))
img_min = int(np.min(neigh))
if img2[i,j] < img_min:
img2[i,j] = img_min
if img2[i,j] > img_max:
img2[i,j] = img_max
return img2
def filter_nl_avg(img, area, val):
img2 = deepcopy(img)
img2 = img2.astype('float64')
x,y = img.shape
kernel = np.ones((area, area), np.float64)
kernel_sum = np.sum(kernel)
offset = (area- 1) // 2
for i in range(offset, x - offset):
for j in range(offset, y - offset):
img2[i,j] = convolution(img, kernel, offset, i,j) / kernel_sum
img3 = deepcopy(img)
for i in range(0, x):
for j in range(0, y):
sum = 0
scalar = 0
for p in range(0, x):
for q in range(0, y):
w_val = math.exp(-(img2[i,j]-img2[p,q])**2)
scalar += w_val
sum += img[p,q]*w_val
img3[i,j] = sum / scalar
return img3
def add_filter(img_filt, fun, area, val):
img2 = fun(img_filt, area, val)
fix, ax = plt.subplots(1,2)
fix.set_size_inches(20,10)
ax[0].imshow(img_filt, cmap='gray')
ax[1].imshow(img2, cmap='gray')
plt.show()
print("RMSE przed filtrem: ", rmse(img,img_filt), " RMSE po filtrze: ", rmse(img,img2))
print("SNR przed filtrem: ",decybel(snr(img, img_filt)), " SNR po filtrze: ", decybel(snr(img,img2)))
print("PSNR przed filtrem: ",decybel(psnr(img, img_filt)), " PSNR po filtrze: ", decybel(psnr(img,img2)))
add_filter(img_salt_pepper_1, filter_avg, 3, 1)
add_filter(img_salt_pepper_1, filter_avg, 5, 1)
add_filter(img_salt_pepper_1, filter_avg, 7, 1)
add_filter(img_gauss_10, filter_avg, 3, 1)
add_filter(img_gauss_10, filter_avg, 5, 1)
add_filter(img_gauss_10, filter_avg, 7, 1)
add_filter(img_salt_papper_gauss, filter_avg, 5, 1)
add_filter(img_salt_pepper_1, filter_gauss, 5, 0.8)
add_filter(img_salt_pepper_5, filter_gauss, 7, 1)
add_filter(img_salt_pepper_5, filter_gauss, 9, 1.5)
add_filter(img_gauss_10, filter_gauss, 7, 1)
add_filter(img_gauss_10, filter_gauss, 9, 1.5)
add_filter(img_gauss_20, filter_gauss, 11, 1.7)
add_filter(img_gauss_50, filter_gauss, 11, 1.7)
add_filter(img_salt_papper_gauss, filter_gauss, 7, 1)
add_filter(img_salt_pepper_5, filter_median, 3, 1)
add_filter(img_salt_pepper_10, filter_median, 3, 1)
add_filter(img_salt_pepper_20, filter_median, 3, 1)
add_filter(img_salt_pepper_20, filter_median, 5, 1)
add_filter(img_salt_pepper_40, filter_median, 5, 1)
add_filter(img_salt_pepper_40, filter_median, 9, 1)
add_filter(img_gauss_5, filter_median, 5, 1)
add_filter(img_gauss_5, filter_median, 7, 1)
add_filter(img_salt_pepper_01, filter_conservative_smoothing, 3, 1)
add_filter(img_salt_pepper_5, filter_conservative_smoothing, 3, 1)
implementacja fft 1d i 2d od zera
def fast(x, inv=0): #impl rekurencyjna, drugi arg odpowiada za odwracanie
n = x.shape[0]
nums = np.asarray(x)
if n == 1: return x
left = fast(x[::2], inv) #wywolanie reku dla parzystych i niepazystych
right = fast(x[1::2], inv)
if inv==0:
inverted_root = np.exp(2 * np.pi * 1j / n)
if inv==1:
inverted_root = np.exp(2 * np.pi * -1j / n)
root = 1
res = np.zeros((n),dtype=complex)
for i in range(0, int(n/2)):
res[i] = left[i] + root * right[i]
res[int(i + n/2)] = left[i] - root * right[i]
root = root * inverted_root
if inv==1:
res[i]/=2
res[int(i+n/2)]/=2
return res
def fft2d(x, inv=0): #fast fourier transform 2D
sh = x.shape
z = np.zeros((sh[0],sh[1]), dtype=complex)
if inv==0:
for i in range(sh[0]):
z[i,:] = fast(x[i,:],0)
for i in range(sh[1]):
z[:,i] = fast(z[:,i],0)
if inv==1:
for i in range(sh[0]):
z[i,:] = fast(x[i,:],1)
for i in range(sh[1]):
z[:,i] = fast(z[:,i],1)
return z
def shift(x): #przeksztlacenie S z prezentacji
n = x.shape[0]; m = x.shape[1]
y = np.copy(x)
y[int(n/2):, int(m/2):] = x[:int(n/2), :int(m/2)]
y[:int(n/2), :int(m/2)] = x[int(n/2):, int(m/2):]
y[int(n/2):, :int(m/2)] = x[:int(n/2), int(m/2):]
y[:int(n/2), int(m/2):] = x[int(n/2):, :int(m/2)]
return y
Lena w częstotliwości
im_lena = Image.open('lena512.bmp').convert('L') # grayscale
lena = np.array(im_lena)
plt.figure(figsize=(8,8))
plt.imshow(lena, cmap='gray')
plt.show()
freq = fft2d(lena)
freq = shift(freq)
plt.figure(figsize=(8,8))
plt.imshow( (20*np.log10( 0.1 + freq)).astype(int))
plt.show()
freq = fft2d(lena)
plt.figure(figsize=(8,8))
plt.imshow( (20*np.log10( 0.1 + freq)).astype(int))
plt.show()
filtry w dziedzinie częstotliwości
def ideal_low(x,r):
r = r*r #r = odleglosc od srodka
M = x.shape[0]; N = x.shape[1]
pt = np.sum(np.square(np.absolute(x)))
p = 0
for i in range(M):
dist_i = (i-M/2)*(i-M/2)
for j in range(N):
if dist_i + (j-N/2)*(j-N/2) > r:
x[i,j]=0
else:
p += np.square(np.abs(x[i,j]))
print('zachowany procent mocy = ', p/pt*100)
return x
def ideal_ell(x,r): #ellptical cutoff
r = r*r #r = odleglosc od srodka
M = x.shape[0]; N = x.shape[1]
coeff = M/N
for i in range(M):
dist_i = (i-M/2)*(i-M/2)
for j in range(N):
if dist_i + (j-N/2)*(j-N/2)*coeff*coeff > r:
x[i,j]=0
return x
def butterworth(x,r,n):
M = x.shape[0]; N = x.shape[1]
pt = np.sum(np.square(np.absolute(x)))
p = 0
for i in range(M):
dist_i = (i-M/2)*(i-M/2)
for j in range(N):
dist_j = (j-N/2)*(j-N/2)
par = 1/(1 + ((np.sqrt(dist_i+dist_j))/r)**(2*n))
x[i,j] *= par
p += np.square(np.abs(x[i,j]))
print('zachowany procent mocy = ', p/pt*100)
return x
def gaussian_low(x,r):
r = r*r #r = odleglosc od srodka
M = x.shape[0]; N = x.shape[1]
pt = np.sum(np.square(np.absolute(x)))
p = 0
for i in range(M):
dist_i = (i-M/2)*(i-M/2)
for j in range(N):
dist_j = (j-N/2)*(j-N/2)
par = np.exp(-(dist_i + dist_j)/(2*r))
x[i,j]*= par
p += np.square(np.abs(x[i,j]))
print('zachowany procent mocy = ', p/pt*100)
return x
def filtering_choice(image, method, r=0, n=0):
freq = fft2d(image,0)
freq = shift(freq)
if method=='low': freq = ideal_low(freq,r)
elif method=='butter': freq = butterworth(freq,r,n)
elif method=='gauss': freq = gaussian_low(freq,r)
else:
print("invalid method name")
return
freq = shift(freq)
rev = fft2d(freq,1).real
plt.figure(figsize=(15,10))
plt.imshow(rev, cmap='gray')
plt.show()
return rev
def f1(t):
return 1/(1+(t/50)**2)
def f2(t):
return 1/(1+(t/50)**4)
def f3(t):
return 1/(1+(t/50)**6)
def f4(t):
return 1/(1+(t/50)**8)
t1 = np.arange(0,200,0.5)
plt.figure(figsize=(8,8))
plt.xlim(0,200)
plt.ylim(0,1)
plt.plot(t1, f1(t1), 'b-', label='n=1', linewidth=2)
plt.plot(t1, f2(t1), 'r-', label='n=2', linewidth=2)
plt.plot(t1, f3(t1), 'g-', label='n=3', linewidth=2)
plt.plot(t1, f4(t1), 'm-', label='n=4', linewidth=2)
plt.title('Przekrój filtru Butterwortha dla różnych n')
plt.legend(prop={'size': 16})
plt.savefig('butter.pdf')
def f1(x):
return np.exp(-x*x/(2*10*10))
def f2(x):
return np.exp(-x*x/(2*20*20))
def f3(x):
return np.exp(-x*x/(2*40*40))
def f4(x):
return np.exp(-x*x/(2*100*100))
t1 = np.arange(0,300,0.5)
plt.figure(figsize=(8,8))
plt.xlim(0,300)
plt.ylim(0,1)
plt.plot(t1, f1(t1), 'b-', label=r'$D_{0}=10$', linewidth=2)
plt.plot(t1, f2(t1), 'r-', label=r'$D_{0}=20$', linewidth=2)
plt.plot(t1, f3(t1), 'g-', label=r'$D_{0}=40$', linewidth=2)
plt.plot(t1, f4(t1), 'm-', label=r'$D_{0}=100$', linewidth=2)
plt.title(r'Przekrój filtru Gaussowskiego dla różnych $D_{0}$')
plt.legend(prop={'size': 16})
#plt.show()
plt.savefig('gauss.pdf')
plt.figure(figsize=(8,8))
plt.xlim(0,200)
plt.ylim(-0.1,1.1)
plt.hlines(1,0,100, colors='b', linewidth=3)
plt.hlines(0,100,200, colors='b',linewidth=3)
plt.vlines(100,0,1, colors='b',linewidth=3)
plt.title(r'Przekrój filtra idealnego dla $D_{0}=100$')
plt.savefig('ideal.pdf')
from scipy.stats import multivariate_normal
from mpl_toolkits.mplot3d import Axes3D
mu_x = 0
variance_x = 8
mu_y = 0
variance_y = 8
x = np.linspace(-10,10,500)
y = np.linspace(-10,10,500)
X, Y = np.meshgrid(x,y)
pos = np.empty(X.shape + (2,))
pos[:, :, 0] = X; pos[:, :, 1] = Y
rv = multivariate_normal([mu_x, mu_y], [[variance_x, 0], [0, variance_y]])
fig = plt.figure(figsize=(10,10))
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, rv.pdf(pos)/2,cmap='inferno',linewidth=0)
plt.savefig('gauss2d.pdf')
plt.show()
moon = (Image.open('moonlanding.png'))
moon = moon.crop((50,200,562,456))
moon = np.array(moon)
moon.shape
lena_salt = add_salt_paper(lena,5)
plt.figure(figsize=(15,10))
plt.imshow(lena_salt, cmap='gray')
filtering_choice(lena_salt, 'butter', 70, 2)
l50 = (Image.open('l50.png')).convert('L')
l50 = l50.crop((30,30,542,542))
l50 = np.asarray(l50)
plt.figure(figsize=(15,10))
plt.imshow(l50, cmap='gray')
filtering_choice(l50, 'low', 30)
filtering_choice(l50, 'low', 40)
filtering_choice(l50, 'low', 50)
filtering_choice(l50, 'butter', 35, 2)
filtering_choice(l50, 'butter', 45, 2)
filtering_choice(l50, 'butter', 60, 2)
filtering_choice(l50, 'gauss', 40)
filtering_choice(l50, 'gauss', 50)
filtering_choice(l50, 'gauss', 60)
zdjęcie księżycowe
plt.figure(figsize=(15,10))
plt.imshow(moon, cmap='gray')
#plt.show()
plt.savefig('moon_org.png')
for i in [120,80,50,30,20,10]:
filtering_choice(moon, 'gauss', i)
filtering_choice(moon, 'low', 40)
filtering_choice(moon, 'low', 50)
filtering_choice(moon, 'low', 60)
filtering_choice(moon, 'butter', 20,1)
filtering_choice(moon, 'butter', 27,2)
filtering_choice(moon, 'butter', 38,4)
filtering_choice(moon, 'gauss', 24)
porównanie na takim samym R
filtering_choice(moon, 'low', 30)
filtering_choice(moon, 'butter', 30,2)
filtering_choice(moon, 'gauss', 30)
wyświetlenie częstotliwości - porównanie
freq = fft2d(moon,0)
freq = shift(freq)
psd2D = np.abs( freq )**2 #moc spektrum
plt.figure(figsize=(15,10))
plt.imshow( np.log10( psd2D ))
plt.show()
ideal_low(freq,60)
psd2D = np.abs( freq )**2 #moc spektrum
plt.figure(figsize=(15,10))
plt.imshow( np.log10( psd2D ))
plt.show()
freq = fft2d(moon,0)
freq = shift(freq)
butterworth(freq,60,2)
psd2D = np.abs( freq )**2 #moc spektrum
plt.figure(figsize=(15,10))
plt.imshow( np.log10( psd2D ))
plt.show()
freq = fft2d(moon,0)
freq = shift(freq)
gaussian_low(freq,40)
psd2D = np.abs( freq )**2 #moc spektrum
plt.figure(figsize=(15,10))
plt.imshow( np.log10( psd2D ))
plt.show()
zjawisko ringingu
im33 = (Image.open('s44.png').convert('L')) # conert to greyscale
im33 = im33.crop((0,0,512,256))
im33 = np.array(im33)
print(im33.shape)
plt.imshow(im33, cmap='gray')
plt.show()
filtering_choice(im33, 'low',20)
filtering_choice(im33, 'gauss',20)